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Abstract. Surface BRDF influences not only radiance just above the surface, but that 
emerging from the top of the atmosphere (TOA). In this study we propose a new, fast and 
accurate, algorithm CASBIR (correction for anisotropic surface bidirectional reflection) 
to account for such influences on radiance measured above TOA. This new algorithm is 
based on 4-stream theory that separates the radiation field into direct and diffuse 
components in both upwelling and downwelling directions. This is important because the 
direct component accounts for a substantial portion of incident radiation under a clear 
sky, and the BRDF effect is strongest in the reflection of the direct radiation reaching the 
surface. The model is validated by comparison with a full-scale, vector radiation transfer 
model for the atmosphere-surface system (Ahmad and Fraser, 1982). The result 
demonstrates that CASBIR performs very well (with overall relative difference of less 
than 1%) for all solar and viewing zenith and azimuth angles considered in wavelengths 
from UV to near-IR over 3 typical, but very different surface types. Application of this 
algorithm includes both accounting for non-Lambertian surface scattering on the 
emergent radiation above TOA and a potential approach for surface BRDF retrieval from 
satellite measured radiance. 


1. Introduction 

The radiance observed by satellites at the top of the atmosphere is the backscattered 
portion of incident solar radiation by the Earth's ground-atmosphere system. This 
radiance can be separated into two components: one from purely atmosphere 
backscatterina (path reflection), and the other from the radiance reflected by the 
underlying surface and transmitted through the atmosphere toward the satellite or in- 
atmosphere instrument. Traditional atmosphere radiative transfer (RT) models usually 
assume isotropic scattering of the lower boundary (Lambertian surface) when calculating 
the contribution from the underlying surface (see Dave. 1964, for example). 

However, natural surfaces are usually non-Lambertian. i.e.. they scatter light 
anisotropically. They all exhibit some degree of bidirectional reflection properties. The 
most common are specular reflection that occurs for forward scattering when incident 
angle is equal to reflection angle (e.g., sunglint for water surfaces), and hotspot 
phenomenon where the strongest backscattering occurs when the source light is exactly 



behind the viewer tor porous media such as vegetated surfaces (Qin and Goei 1995) The 
change surtace reflectances with bolh s „lar and viewing directions , s often refeej » 
as the surface BRDF (bidirectional reflectance distribution function) property This 
property has been considered and investigated by simulation of radiation propagation in 

— covers Jde^ 

Surface BRDF affects the emergent radiation from the top of the atmosphere as well 

(Coulson‘“ T 1 966)' and' TlTd b °‘ h “ Ra - vlei - ah atmosphere 

1 r IS" 966) and a turbld atrnos phere (Keopke and Kriebel 1978) found 
significant differences in radiance a, TOA over natural surfaces and their £amben-modet 
quivalents even though their albedos were equal. Fitch (1981) did a similar studv for 
three types of natural soil surfaces based on laboratory m eSurem”« of 
bidirectional reflectance combined with an atmosphere-surface model He considered 
polarization and multiple scattering between the two media (the atmosphere and the 
underlying surface) up to five orders of scattering in the vector radiative transfer model 
fonad and Fraser (1982) developed a full-scale, vector RT model for the atmosphere- 
ocean system, in which the anisotropy of scattering from a rough ocean is incorporated. 

surface ' ^anre" the ^erlying non-Lambertian 

surtace lanre et al. (198 j) specified the radiation emerging from the top of the 

atmosphere (TOA) in five pans: incident direct/diffuse radiation, directly/diffusely 

transmitted through the atmosphere after surface reflection, plus a term for multiple 

the eri n 8 ? etVVeen the SUrfaCe and the atmos phere. In their final formulation, however 
he multiple-scattering contribution is ignored. The 6S model [Vermote et al. 19971 

a e a term to the above formulation trying to approximate multiple scattering 
contributions. It was not introduced in a physically consistent manner (the authors called 
the arbitrary addition an approximate term in the paper) because there were no 
physically -based derivations to justify this addition. Further, the direct-to-diffuse (or 
diffuse-to-directional) conversion of radiation caused by surface reflection was 
improperly simulated. In the DISORT code, Stamnes et al. (1988) assumed the surface 
BRDF to be a function only of the phase angle, so that Legendre polynomials can be used 
to represent non-Lambertian surface reflection. Since it only represents a single anafo 
hib assumption is not valid in general. Even for a single angle, the ability of Legendre 
functions to represent the surface BRDF is quite limited. 

To remove formulation uncertainties in the application of such algorithms as applied 
o general situations, it is desirable to base the derivation on a consistent consideration of 
the physics. The new algorithm presented below will improve accuracy when accounting 
or the surface BRDF, as well as providing fast computational implementation^ 
Comparisons with a full radiative transfer treatment of the BRDF effect (using a modified 
version of the Guass-Seidel approach based on Ahmad and Fraser. 1982 ) validate the 

algorithm s accuracy and show the need for computational speed in a hfohlv accurate 
approximation. - - 


Radiation incident on a surface is partially specular (collimated) and partially dittu.se. 
The specular component results from transmission of the direct solar beam through the 
atmosphere, while the diffuse component results from atmosphere scattering ot the 
incident solar beam and downward scattering of upwelling radiation after reflection from 
the surface. The BRDF effect is strongest in the reflection of the specular component, 
which can be more than 80% of total incident solar radiation (in visible and near-infrared 
resions) under a clear sky. In order to capture this BRDF effect, the direct (or collimated) 
component has to be treated separately from the diffuse component when modeling its 
interaction with the surface. Also, to completely take account of surface BRDF effects on 
measurements above TOA, one has to consider the conversion between direct or 
directional radiation and diffuse radiation from surface reflection, as well as multiple 
scattering between the atmosphere and the surface. These goals can be achieved by 
applying well-known four-stream theory (Hapke, 1981; Li et al., 1996; Qin and Liang, 
2000) to simulate radiation interaction in the boundary between two media. 

In this study, we will express radiation flux in terms of diffuse (in upwelling and 
downwelling directions) and direct or directional (collimated in a specified direction) 
components. This division will ensure that a correct surface reflection coefficient is 
applied for a given set of incoming and outgoing radiation from surface reflection (see 
section 2.2.1 for details). When dealing with multiple scattering between two media, the 
conversion and multiple interactions among all types of radiation (downwelling and 
upwelling direct and diffuse) are fully taken into account with proper reflection 
coefficients applied, including all orders of multiple scattering. Within the limitations of a 
4-stream approximation, this will overcome most of the weaknesses in previous 
atmosphere-surface models incorporating surface BRDF characteristics in the lower 
boundary. The new algorithm will produce the most complete high-speed approximation 
of multiple scattering between the surface and the atmosphere to date. 

In the following section, we will describe the development of this new algorithm: 
Correction for Anisotropic Surface Bidirectional Reflection (CASBIR), including 
determination of the various coefficients. Then, we will compare the results with that 
from a precise, vector-based atmospheric radiative transfer (R.T) model for a pure Raleigh 
atmosphere [Ahmad and Fraser, 1982]. Some concluding remarks and application issues 
will be discussed in the last two sections. 


2. Algorithm Development 


For a Lambertian surface with surface reflectance r 5 . the reflectance above the 
atmosphere can be expressed as follows 


tf,(v) = D 


r (/) • r ■T'(v) 

1 -Vr 


( 1 ) 


where /. v are illumination (solar) and viewing directions, r o is the path scattering- 
reflectance of the atmosphere. 71 is the total transmittance from the top ot the atmosphere 
to the izround along the path of the incoming solar beam, and T is the total transmittance 
from the ground to the top of the atmosphere in the view direction ot the satellite. St, is 



££?£ the “ pw ; rd ditruse . tlllx ^scattered from the atmosphere to the Earth's 
, ' . l r a non *. anlb °nian surface. surtaee reflectance chances with both illumination 
and viewing directions. This, combined with the anisotropic diffuse irradian IS “m 
upon the surtaee. means Eq.fl, could result in considerable errors in rndian" 

below 8 ° n ' hm “ aCC °“ m f ° r no "- Lamb «i a n surface reflection is derived 

Consider an atmosphere bounded below by a non-Lambertian reflectin« surface with 
bidirectional reflectance rO', e). We divide the radiation field m the medium into mo 
components (fluxes): diffuse (£) and direct (or directional) (F). Symbols "T” and “l- 
s and tor upwellmg and downwelling components with superscript (or subscript) for 
q antities at the upper (or lower) boundary of the atmosphere (see Fig 1) We alsodefine 
a as the path reflectance, , as the path transmittance, and r as the 

h T tW ° subscnpt s y mbols ^ “d" (direct or directional) or 4" [diffuse or 
misphenc (i.e., the average over the hemisphere)], to indicate photon status before and 

m S Ct, °H 1“’ ^ ^ f ° Ur COmbinatio - of these two s" W’ 

’ . and hh v^ith the first symbol indicating the initial status of nhotnnc 
(incoming) and the second one for the resulting photon status after interaction (outgoing) 

beinp"fTfFi P S | Subs <f lpt lndlcates th ^ incoming photons from a specific direction 
being diffusely scattered into the whole hemisphere, either by backward atmospheric path 

£u„d™y reflecfionlLt " f0r " ard ^ ^ W - 

• the ?°T “1 Fi i U at th t t0p 0f the atmos P here > die only incident radiation flux F l (i) 
is the direct solar beam with zenith angle 6, and azimuth angle <p„ The radiation scattered 

ret”^ ° f VieW (F0V) ^ ^ V!eW “ <*• ^ m is the sutTof 

F ‘ ( y ) v ) ' F l (/) + t dd (v) • F. (v) + t M (v) • E r , (2) 

where Fj(v) is the radiation flux in direction v, and Er is the upwelling diffuse flux 
leacmg the surface. a dd is the purely atmospheric backscattering coefficient [also called 
path reflectance, equal to r 0 in Eq.(l)], t dd is the direct part of T\ i e ■ atmospheric 
transmittance for collimated radiation in the satellite viewing direction, and t hd is the 

radial' 6 ^ ° r’ tbe efpIcien cy of atmospheric scattering of upward diffuse 

radiation coming from the surface toward the satellite in the viewing direction (* v . *.) at 

TOA. W e will discuss these and other coefficients in details in Section 2.2. 

2.1. Calculation of F-(v) and E- 

Based on the above definitions, the equations to compute F- and E - can be expressed as 

F- ( v) = r M (/. v) • F (/) + r, a (v) E 

£ * = r , i A‘)-Fji) + r, i ,-E. (3a) 

Tb " f b ° h Ve eq “ ations 5 c °nsider the conversion between directional and diffuse radiation 
caused by surface reflection. The solution of the above equation set can be obtained bv 
first evaluating the 1 ’'-order scattering components of E and F, then the "--order 
components, the a -order components, and so on. By considering multiple bounces of 
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photons between two layers/media, Li et al. (1996) deduced a set ot closed form 
expressions including all orders ot multiple bounces tor the above quantities. Based on 
our definitions, the solution to the above equations tor the lower boundary ot the 
atmosphere-surface system can be written as follows: 


/> (v) = r Jd (/, v) ■ F i (/) + /*, (v) ■ [r Jh (/) -F*(i) + <J 


hh 


•£.J 


£r = 


r Jh {i)-FM) + r, 




0. 


1 - r., ■ a ■„ 


(3b) 


F l W = t„W-F l V)- 

replacing F^ and E j in Eq.(2) with the above expressions in 3a and 3b, and making 
some mathematical manipulations, we finally get the reflectance at TOA — R a , defined as 
the ratio of F\v)/F^(i), in a form of matrix, as 


T(/) • R(/, v) • T(v) - t Jd (/) • t dd (v) ■ |R(i, v)| • <7 hh 

KV, v) = .o M O', v) + 1 


(4a) 


where matrices T(i), T(v) and R are defined as 


T(0 = L(0 f*(0l T(v) 


'*00 

M V )J 


R(»>) 


r dd (hv) r lh (i ) 
hjiy) r hh 


(4b) 


One can see Eq.(4a) is similar to Eq.(l) in form, with 

Ti(i)=tjd(i) + tdh(i), V(y)=t da {v)+ tUv), and 

F 0 = C7 ddi Sb~F5hh- F s F),h- 


(4c) 


For non-Lambertian surfaces, single variables such as T\ T l and r s in Eq.(l) are replaced 
by the corresponding matrices. Also, there is an extra term in Eq.(4a), w'hich is a function 
of the determinant |R|, calculated as 


|R| r dd'F hh~Fdh'fhd* (4d) 

Since |R| could be positive or negative, depending on solar and viewing directions, and 
the degree of non-Lambertian reflection from the surface, the contribution trom non- 
Lambertian surface reflection could be more or less than that from its Lambert 
equivalent. For a Lambertian surface, Eq.(4a) reduces to Eq.(l). because the four 
components in R are equal to f m „ so that |R|=0. 


2.2. Estimation of various coefficients 

There are two types of coefficients involved in CASBIR [Eq.(4a)]: atmosphere 
related [T(/). T(v). each has two components, and a. which has tour components, but 
only two (ojj, a hh ) are used here] and surface related (R. which has four components). 
Fortunately, as shown below, coefficients in one group are independent ot those in the 
other group, and can be determined separately. For example, atmospheric-property 



related coefficients could be calculated by using any atmospheric RT model based on a 
Lambertian surface assumption. Similarly, components ofR can be 

^^B^Ff eCtar : Ce dlStnbUtIOnS - regardleSS ° f at m osphe“ dl t!r 

ause surface BRDF is an intrinsic property of the surface, independent of atmosDheric 
co„d,,,o„, in the following, we win describe delegation o' coeffiden^Tn' bZ 

2.2. 1. Boundary reflectance coefficients 

Surface BRDF is the physical quantity to characterize surface reflection It is 
determined by surface structures and optical properties (material reflectivity and 

is def hT ’ and V r aneS Wlth lllumination and observation directions. However BRDF 
defined for an infinitesimal solid angle; it can be modeled, but virtually cannot be 
measured. In practice, it is often replaced by a measurable alternative Bl^ 
(bidirectional reflectance factor), quantified as the ratio between radiance reflected from a 

~ “ e m 3 PerfeCt Lambertm refleC,0r ' ■“ 

All I components of R are functions of surface BRF and can be calculated 
strmghtfomardly once the surface BRF is determined. In the following, we will discuss 
the definition and determinations of components of R. We leave the discussion on 
practical methods to determine surface BRF to Section 4. 

r d d - bidirectional reflectance r dd can be defined in terms of the BRF, r(i, v), 

r da{i v) =/•(/', v), 

r{i, v)=7i/t(/, v)/j jijF’i(i') 


where 


(5a) 

(5b) 


is the surface BRF under direct solar beam (without any skylight). F l{ i) is the direct solar 
flux incident upon the surface [see Eq.(3b)] in direction / with p, =cos(0,), and /.(/ v) is 
e reflected radiance of F^i) from the surface in direction v. As indicated in section 4 
r(i,s) can be determined by modeling or measurements for a given surface type. 

r . dh ~ di t rectional -hemispheric reflectance r dh specifies the fraction of direct radiation 
incident upon a surface that is diffusely reflected toward the upper hemisphere Pf) 
Mathematically it is defined as i pne 1 h 


r j„ = 


M,F t (/') 

where c is the scattering direction. Replacing / r (/,v) with r(i.v) in Eq.(5b) yields 

r uh = r n (0- r „ d ) = — d(fl: | r{ i. g ) u, du , . 

r;,(i) is the hemispheric reflectance for a specific solar direction. 


(6a) 


(6b) 
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r:,j — hemispheric -directional reflectance r;,j is defined as the traction ot downwelling 
diffuse radiation reflected toward the specific direction v by the surface, i.e.. 


n,J v ) 


| L r (d .v)p.dd. 


(7a) 


where 27t‘ represents the lower hemisphere, C is the incident direction of diffuse light 
upon the surface, and L\ is the upward reflected radiance of E which equals r(i,v)E Ju. 
Therefore, we finally have 


r w (v) = r*(v), r.{v) = - V dtp. [r{fv)p ; dp ; . (7b) 

r h {v) is the hemispheric reflectance for a specific viewing direction. If the surface 
reflection follows the reciprocity law, i.e., r(i,v)= r(vj), then r/,(i)= r/,(v) and r,j’ n ~ fhd if 


i=v. 


rhh _ hemispheric-hemispheric reflectance r^h is also called bi-hemispheric reflectance. 
By definition, it is the double integral of the scattered diffuse radiation (I r ) over viewing 
(upper) and illumination (lower) hemispheres divided by the downwelling diffuse 
radiation. That is. 


r hh _ 


L {L-MooMn- k d£2 -- 


(8a) 


where £, are source light and scattering directions, respectively. Similar to calculate r^, 
the above double integration can be simplified and evaluated as 


r hh = albedo = — f " dtp , | r h {f)p 4 dp, . (8b) 

where r h {p) is the same as r h (v), expressed in Eq.(7b). Therefore, r hh is the spherical 
albedo that considers all viewing and illumination directions. 


2.2.2. Atmospheric scattering and transmission coefficients 

Atmospheric path scattering (a) and transmission (T) coefficients are functions of the 
atmospheric optical depth (t a ), single scattering albedo (to), and phase function ( P ) of the 
scattered and absorbers in the atmosphere. Only the direct component of T [see Eq.(4c)] 
has an analytical expression as 

WH- *a)= exp(-T a /p), (9) 

where p=cos(9). 0 is the zenith angle of the light. To estimate other components of a and 
T, one has to utilize an atmospheric RT model, because there are no analytical 
expressions available in general. Most atmospheric RT models can provide these 
coefficients in a form of either numerical solutions or look-up tables (LLTs) for a variety 
of sun-view geometries and aerosol loadings (e.g.. Dave and Gazdag, 1970). The LET 
approach is computationally more efficient, especially tor practical use in an algorithm 
for satellite fields of view in a global data set. It has been used for decades (e.g., Dave et 
al., 1966). However, under some special conditions (e.g., tor a pure Rayleigh 
atmosphere), one may be able to obtain an analytical approximation for these coefficients 
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mim™° tL | t,n f Tjnrc- . ,99 “^ These analytical expressions become closer to exact 
uimcncal solutions tor wavelengths much larger than the particle size of >nseous 
constituents, because the multiple scattering contribution becomes verv small Since the 

diff^p n g C ° t ‘t nt fpatH refleCtance) and (backward scattering of the upward 
tfuse flux from the surface) are the same as r 0 and in Eq.(l). their determination will 

not be discussed here. In the following, we will only discuss the determination of the 
diffuse components of atmospheric transmission function T. 

%^ l T ti0nal ' h TAt herk path lransmltance r dh defines the fraction of downward 
ltfuse flux generated by atmospheric scattering as direct solar beam passes through the 
atmosphere It is also called the atmospheric diffuse 
Mathematically, it can be expressed as follows 

l dh\ l ) - — 7 , 

M,F l ( 0 ( 10 ) 

where Z I is the downwelling diffuse radiance reaching the surface due to atmospheric 

scattering of the incident direct solar radiation from direction / into direction <f. Generally 

there is no analytical expressions for Z, because of the multiple scattering Tn the 

atmosphere. However, numerical results for Zj are available for a given atmosphere tvpe 

nl!Tr COn th entl °? at " 10s P heric RT mod ^s. Therefore, is usuafly provided for aiven 
solar zenith angles and atmospheric conditions (e.g., a given set of x a , co and P). 

tM - hemispheric-directional path transmttance t hd is defined as the fraction of upward 
fuse flux scattered by atmospheric constituents (molecules and aerosols) toward the 
satellite m direction v. Similar to t d h, thd can be estimated from 


Zk/( V ) = 




( 11 ) 


where Z T is the upwelling diffrise radiance at the top of the atmosphere scattered toward 
viewing direction v by atmospheric constituents. With the same reasons for L Z r can 

evlated n . UmeriCally C ° mPUted m ° St atm0S P heric RT models ^ from which t hd can be 


3. Validation 

ro validate CASBIR. we compare the modeled reflectance at TOA with calculations 
from a full-scale vector atmospheric RT model (Ahmad and Fraser. 1982). here called 
model. The reason for using this atmospheric model as the standard for comparison 
because of its ability to directly incorporate arbitrary surface BRF distributions into the 

T) t eded a m cVsBI^ f ^ Ae » Parameters (a and 

I) needed in CASBIR tor validation purpose. A Ravleiah atmosphere in five 

wavelengths (388 nm. 443 nm, 551 nm. 645 nm, and 870 nm. matchina 5 of 10 Triana 

channels) covering L V to near-IR over three diverse types of surface (desert Grassland 

and torest) are considered in the comparison. The VRT model has been successive 

compared with both DISORT and the Dave vector code for Lambertian surfaces 
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To obtain surface BRF distributions, an elaborate 3-D scene BRDb model (Qin et al., 
1 998; Qin and Gerstl, 1998; 2000) is employed to first generate the 3-D structures of a 
given surface type and then compute the complete BRF distributions. The input optical 
properties (reflectance and transmittance of vegetation elements and the soil background) 
and structural parameters are taken from field measurements (Privette et al., 2000; 
Waiter-Shea et al., 1992; Hall et al.. 1992). Our choice of modeled surface BRFs with the 

3- D scene model rather than using field BRF measurements is based on the high angular 
resolution of the simulated BRF data. This enables us to complete the comparison over 
the whole hemisphere without having to do any interpolations for surface BRF. 

Specifically, for each surface type in each wavelength, we produced a look-up table 
for surface BRF over 16 view zenith angles (VZAs, 6° step) and 16 relative azimuth 
angles (12° step) for each solar zenith angle (SZA), with a total of 16 LUTs generated for 
16 SZAs. Figure 2 shows surface BRF distributions for 551 nm at a SZA of 30° for 3 
surface types. Note that the desert exhibits the largest contrast in reflectance between 
forward and backward directions because its rolling-hill structure (which reflects the 
topography of natural deserts) produces a significant amount of shadowing near the 
forward scattering directions. Forest scenes have the sharpest hotspot reflection peak 
because of the highly heterogeneous structure (large gaps exist in the canopy), and the 
grassland has the least variations in reflectance compared to the above two. Both 
grassland and forest have flat soil backgrounds. 

Based on the LUTs for surface BRF, the VRT model calculates reflectance at TOA 
for the Rayleigh atmosphere overlying each non-Lambertian ground surface. It also 
produces other parameters needed for CASBIR, such as path scattering (errand G/,/,) and 
transmission (t dh and t hd ) coefficients. Figure 3 plots the corresponding BRFs at TOA for 
each surface BRF distribution shown in Fig. 2. Comparison between Figure 2 and 3 for 
all five wavelengths (not shown here) indicates that BRF effects increase with longer 
wavelengths, because of the smaller atmospheric optical depth and smaller path- 
scattering contribution. 

Comprehensive comparisons of results are presented in Figs 4-6 and Table 1. Figures 

4- 6 examine the performance of CASBIR in the solar principal plane (view azimuth in 0- 
180° transect) for the three types of surfaces. We also include the result from the Lambert 
equivalents for comparison. For each surface type, we plot three wavelengths (388, 551 
and 870 nm) under three solar zenith angles (6, 30. and 60’). CASBIR performs very well 
for all cases, matching the distributions from the VRT model at all points. As a contrast, 
however, the Lambert equivalent [Eq.(l) with r s = f"hh] produces substantial differences 
(except for the ultraviolet band, such as 388 nm). especially in the hotspot region or in the 
visible and near-IR. 

Table 1 summarizes the mean and maximum relative differences between the VRT 
model and CASBIR (or its Lambert equivalent) for all solar and viewing directions 
considered, and for each surface type (a total of 2912 cases after excluding zenith angles 
larger than 80 J ). Numbers in table 1 prove that this simple algorithm is very effective in 
accounting for surface BRDF influences at all angles. The average percentage difference 
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is below 0.5 o in the L \ tor all surface types, and over vegetation cover for all 
vvave <-ngths from I V to near-IR. However, the relative difference goes up to 1% for the 
desert m the visible. The mean percentage difference under Lambertian assumption is 
much higher, with the highest up to 43% for the desert in the near-IR. Generally, the 

nDnu" 0 ^ ' S larg ? L aS SUrface reflectance increases, indicating the enhanced surface 
B , R P F ' nfluence ' Therefore, the surface effect is wavelength dependent, simplv because 
ot 1) changes in surface material reflectivity with wavelength, and 2) the wavelength 
dependence ot atmosphere scattering. The percentage difference also changes with solar 
zemth angle, increasing as SZA decreases (see Figs. 4-6). because the contribution from 
atmospheric path scattering decreases with SZA. Finally, for a given surface type at a 
lxe wavelength, the surface BRDF effect varies with the type of surface formations and 
topography, since highly heterogeneous surfaces usually have strong anisotropic 
scattering, and accordingly, strong surface BRDF. 


4. Discussion 


For visible and near-IR radiances under clear skies, the above studies demonstrate that 
surface BRDF has a considerable influence on radiation emerging from the top of the 
atmosphere. However, evaluation of such effects requires accurate or reliable information 
on spatial distribution of surface bidirectional reflectance. Currently, this information is 
mostly obtained from two sources: field measurements or BRDF model simulations. 
Although field measurements can directly provide the ground truth of surface BRF, it is 
very labor intensive and only available for very limited areas and a few types of land 
cover. Need for surface BRF at a global scale cannot be met with this method. On the 
other hand, BRDF models have the capacity to generate BRF distributions globally. But 
they have their own limitations: inflexible model applicability and difficulties in 
determining the needed model input parameters. 

Most BRDF models (except 3-D models) only work over specific surface types and 
their input parameters are not obtainable globally. Therefore, to solve the above 
problems, the MODIS surface BRDF/albedo retrieval team uses a kernel-driven model (a 
linear sum of pre-specified terms characterizing different scattering modes) combined 
with model inversion technique to obtain the input model parameters to reconstruct 
surface BRDF (Strahler and Muller. 1999). First, it is assumed that the kernel-driven 
model has a universal applicability for all the surfaces the satellite will observe. Then, 
satellite measurements are used to retrieve the model parameters so that the BRF 
distribution can be calculated by running the model with the retrieved parameters 

Even though the kernel driven method seems the only practical choice to make use of 
satellite data, this approach still faces a few challenges. First, although kernel-driven 
models are simple and fast, they are not universal, because land surfaces are very diverse' 
its scattering nature cannot be characterized by a linear sum of two or three simple 
predetermined kernels. Therefore, the model may not have the capacity to capture the 
BRDF ot every type ot surface the satellite observes (particularly for hfrhlv 
heterogeneous scenes or mixture pixels). Second, the input parameters retrieved from 
satellite measurements may not be unique, which could lead to the BRDF pattern 
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reproduced by such models that is different from the true one for the surface viewed in 
the pixel. Third, the satellite data used for surface BRDF'albedo retrieval have to be 
atmospherically corrected. The dilemma here is that the algorithm tor atmospheric-eftect 
correction could not produce a correct result for non-Lambertian surfaces without tirst 
knowing surface BRF distributions. Therefore, further work is still needed in order to use 
satellite data for reliable surface BRDF retrieval. 

Besides accounting for surface effects on radiation emerging from TO A. our new 
algorithm (CASBIR) can also be used to retrieve surface BRDF from satellite 
measurements. Theoretically, this can be done by solving the integral equation [Eq.(4a)] 
with multi-angular satellite measurements for a clear sky or if the atmospheric profile is 
known. But practically, there is no guarantee that a unique solution exists for surface 
BRDF from Eq.(4a). For most satellites there are usually not enough angular samples to 
ensure a reliable retrieval of surface BRDF from satellite measurements in a short period 
during which surface BRDF does not change much. These constraints make the direct 
retrieval approach less attractive and useful. 

How'ever, some new thoughts about surface BRDF retrieval can be gained from the 
result of CASBIR. For example, the similarity of angular patterns of reflectance between 
above-the-surface and above-the-atmosphere in the near-IR (see Fig. 7) indicates that the 
angular distribution pattern of surface BRDF is well retained in satellite observations for 
near IR wavelengths under a clear sky. Therefore, if we can obtain such patterns from 
other sources, we can retrieve surface BRDF from just a single-direction satellite 
observation (e.g., from nadir). Obviously, a dedicated 3-D BRDF model can help to 
provide such BRDF patterns if reliable global or regional land cover maps are available. 
We will detail this approach in other papers in the context of using Triana observations 
combined with measurements from other satellites (such as SeaWiFS, MODIS, MISR, 
etc.) to estimate surface radiance and energy budget. 

5. Conclusions 

A new, fast algorithm to account for non-Lambertian surface scattering on radiation 
emerging from the top of the atmosphere is developed in this study. Rather than treating 
photons in every direction equally and precisely, as in the VRT model, in its interaction 
with the surface, we group radiation into direct and diffuse categories and treat both 
groups separately. The physical basis for this lies in the fact that the BRDF effect is the 
strongest in the reflection of the direct incident radiation, which comprises a substantial 
proportion in the incident radiation under a clear sky. The separation allows us to apply 
four-stream theory to handle the complicated problem of radiation interaction between 
two media (the atmosphere and underlying surface), and to develop a simple, analytical 
algorithm to account for surface BRDF effects on satellite measured radiance. 

The comparison with an accurate, full-scale vector atmospheric radiation transfer 
(VRT) model that is capable of directly incorporating arbitrary surface BRDF as its lower 
boundary condition, demonstrates that the new fast algorithm is very accurate and 
effective. The relative difference is less than 0.5% meanly in the LA' region for all three 


surface types (desert, grassland and forests) or in the spectral region from UV to ne tr-IR 

° r «*“•'«>■ Only for deserts in the visible, the relattve difference eoe up V % 
average due to largest surface reflection. w F ' 


The surface BRDF effect is wavelength dependent. It decreases with increasing 
atmospheric optical thickness and surface reflectance. For a clear skv (a Ravleigh 
atmosphere), the surface influence increases with wavelength due 'to decreasing* 
atmospheric path scattering and increasing surface reflectance. For example in the UV 
reg.on. the contribution of atmospheric molecular scattering dominates, and surface 
reflectance is very small (less than 3% tor a vegetation surface and 7% for desert or bare 
soil surfaces), thus the surface BRDF effect is marginal and can be neglected As the 
wavelength increases, the contribution by atmosphere scattering declines and' surface 
reflectance nses Therefore, in the near-IR, the contribution from surface reflection 
dominates and TOA BRF has almost the same shape as the surface BRF (see the 
Discussion section). This suggests a new approach to retrieve surface BRDF patterns 
from satellite observations in the near-IR. F 
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Tabic 1. The percentage differences of CASBIR and Lambertian model 
for all solar and viewing directions considered (excluding 
zenith angles larger than 80°) over three different surface types 
under a clear sky (Rayleigh atmosphere). 


wavelength 388 443 551 645 8. 0 

Desert (nm) 

CASBIR mean 0.41 0.87 1.22 1.00 0.43 

maximum 1 .00 2.25 4.36 4.58 2.67 

Lamb. mean 2.10 6.73 19.83 

maximum 7.36 18.42 69.73 


Grassland 


CASBIR 

mean 

maximum 

0.06 

0.27 

0.12 

0.62 

0.65 

2.05 

0.13 

0.71 

0.42 

1.17 

Lamb. 

mean 

0.56 


7.45 

8.42 

19.15 


maximum 

4.52 


32.75 

37.90 

66.21 

Forest 

CASBIR 

mean 

0.11 

msm 

0.24 

0.23 

0.28 


maximum 

0.39 

MSM 

0.91 

0.72 

1.47 

Lamb. 

mean 



8.37 

11.86 

13.22 


maximum 



41.07 

54.53 

69.31 


Lamb. 







Figure captions 


F[gure 1 A sketch graph to illustrate radiation interaction in atmospheric 
using the 4-stream scheme (see text for details). 


boundaries 


Figure 2. Simulated surface BRFs at solar zenith angle of 30° in 551 nm over (a) desert 
(b) grassed and (c) forest. The polar coordinate system represents view zenith an* 
• h 0 ( n£ j dir ) at the center of the plot and 90° at the edge. The solar azimuth angle 
increases clockwise with the hot spot direction at 180° and forward scattering direction 


Figure 3. Same as in Figure 2 but for reflectance at the top of the atmosphere. 

FlgUre , 4 . Comparison among TOA reflectances over the desert scene from the VRT 
model (solid line), CASBIR (plus sign) and the Lambertian equivalent (dash line) in 

551 “ d 870 "" (,0P '° b0tt0m) at 


Figure 5. Same as in Figure 4 but over grassland. 
Figure 6. Same as in Figure 4 but over forest scene. 


Figure 7. Similarity of BRF shapes between above-the-surface and above 
at solar zenith angle of 30° in 870 nm over (a) desert, (b) grassland and 


-the-atmosphere 
(c) forest. 
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Up per boundary 




Atmosphere 



Figure 1 . A sketch graph to illustrate radiation interaction in atmospheric boundaries 
using the 4-stream scheme (see text for details). 
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